(N 
O 
O 
(N 

O 

(D 

Q 
o 

(N 

r 

o 
o 

O 
• 1— I 

C/2 



(N 
oo 
O 
(N 

(N 
O 

o 

• 1— I 
Ph. 



X 



Extended Recurrence Plot Analysis and its 
Application to ERP Data 

Norbert Marwan^ * Anja Meinke^ 

^ Nonlinear Dynamics Group, Institute of Physics, 
University of Potsdam, Potsdam 14415, Germany 

^ Institute of Linguistics, 
University of Potsdam, Potsdam 14415, Germany 

September 18, 2002 



Abstract 

We present new measures of complexity and their application to 
event related potential data. The new measures base on structures of 
recurrence plots and makes the identification of chaos-chaos transitions 
possible. The application of these measures to data from single-trials of 
the Oddball experiment can identify laminar states therein. This offers 
a new way of analyzing event-related activity on a single-trial basis. 



1 Introduction 

Neurons are known to be nonlinear devices because they become activated 
when their somatic membrane potential crosses a certain threshold [Kan- 



del et al., 1995]. This nonlinearity is one of the essentials in neural mod- 



elling which leads to the sigmoidal activation functions of neural networks 



[ Amit, 19891. The activity of large formations of neurons is macroscopically 
measurable as the electroencephalogram (EEG) at the human scalp which 
results from a spatial integration of postsynaptic potentials [Nune^ 1981]. 
However, it is an unsolved problem whether the EEG should be treated 
as a time series stemming from a linear or a nonlinear dynamical system. 
Applying nonlinear techniques of data analysis to EEG measurements has 
a long tradition. Most of these efforts have been done by computing the 



correlation dimension of spontaneous EEG [e. g. 


Babloyantz et al., 


1985; 


Rapp et al. , 


1986 


; Gallez and Babloyantz, 1991; 


Lutzenberger et al. , 


1992; 


Pritchard and Duke, 1992 


|. rheiler et al 


[ 


1995 


] applied the technique of 



surrogate data to correlation dimensions of EEG and reported that there is 
no evidence of low-dimensional chaos but of significance for nonlinearity in 
the data. 
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While correlation dimensions are only well defined for stationary time se- 
ries generated by a low-dimensional dynamical system moving around an 
attractor, these measures fail in investigating event-related brain poten- 



tials [ERPs, Sutton et al., 19651 since they are nonstationary by definition. 



Traditionally, ERP waveforms are determined by computing an ensemble 
average over a collection of stimulus time locked EEG trials. This is based 
on the following assumptions: (1) the presentation of stimuli of the same 
kind is followed by the same sequence of processing steps, (2) these pro- 
cessing steps always lead to activation of the same brain structures, (3) 
this activation always elicits the same p attern of elect rophysiological activ- 
ity, which can be measured at the scalp [ Rosier , 198^ and (4) spontaneous 



activity is stationary and ergodic [ beim Graben et al. , 200C]. 



By averaging the data-points time-locked to the stimulus presentation 
(cf Oddball experiment) it is possible to filter out the signal (ERP) of the 
noise (spontaneous activity). In the next step the functional significance 
of a component is assessed. Antecedent conditions of the occurrence of a 
component and variables, which influence its parameters are defined. Now 
the commonalities of these factors are identified. The generalization of all 
empirically found influencing factors leads to a more abstract cognitive the- 
ory of the functional meaning of a event-related potential component and 
makes it usable for the validation of models of cognitive processes. 

The disadvantage of the averaging method is the high number of trials 
needed to reduce the signal-to-noise-ratio [Kutas and Petten, 1994]. This 



is crucial for example for clinical studies, for studies with children and for 
studies, in which repeating a task would influence the performance. So it is 
desirable to find new ways of analyzing event-related activity on a single- 
trial basis. Applying nonlinear methods to electrophysiological data could 
be one way of dealing with this problem. 



To compute dimensions of ERPs, [Molnar et al.| [ 119950 used the pointwise 
dimensions and reported a drop of the pointwise dimension as a function 
of time corresponding to the P300 component observed in the Oddball ex- 
periment. Recently, concepts of information theory have bee n introduced to 



analy se ERPs. On one hand this is the wavelet entropy of iQuiroga et al 



P2001J ] and on the other hand symbolic dynamics of EEG and ER P [beim 
Graben et al, ^OOOt [Frisch et al.] , ^002| ; [Steueij , ^002^ |Schack| , |2002| ]. 



A further promising approach is the recurrence quantification analysis (RQA), 
which is based on the q uantification of the diagonal or iented structures in 
recurrence plots [RPs, Webber Jr. and Zbilut , 1994 ; Zbilut and Webber 
Jr., 1992 ]. The RQA wa s broadly appli e d in a wide field of the analysis of 
physiological data [e. g. pasdagli , 1997; Faure and Korn, 1998; Thomasson 



et al., 2001; Marwan et al., 2002]. The important advantage of methods 



based on the quantification of RPs is that the required data length can be 
relatively short. However, the measures of the classical RQA are only able 
to recognize transitions between periods and chaos and vice versa [Trulla 



et al., 1996]. In this work, we will use recently introduced additional mea- 



sures based on RPs in order to find chaos-chaos transitions in physiological 
data. These new measures use the vertical structures in the RP and are 



able to identify laminar states [Marwan et al., 2002]. 
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In the first section we will give a short introduction into RPs and their 
quantification analysis. In the next section we will introduce the new mea- 
sures and finally we will apply them to event related potential data gained 
from the Oddball experiment. 



2 Recurrence Plots and Their Quantification 



The method of recurrence plots (RP) was introduced to visualize the time 
dependent behavior of the dyn amics of system s , whic h can be pictured as a 
trajectory in the phase space [Eckmann et al., 1987]. It represents the re- 



currence of the m-dimensional phase space trajectory e TZ'" (i = 1, . . . , N, 
time discrete) to a certain state. The main step of this visualization is the 
calculation of the N x A^-matrix 



R 



,J = 1. 



■ N, 



(1) 



where is a state dependent cut-off distance, || • || is the norm of vectors, 8 
is the Heaviside function and N is the number of states. The phase space 
vectors for one-dimensional time series Ui from observations can be recon- 
structed with the Taken's time delay method = {ui,u.i+ r, • • • , Ui+(m-i) r) 



with dimension m and delay r [Kantz and Schreiber, 1997]. The recurrence 



plot exhibits characteristic large-scale and small-scale patterns which are 



caused by typical dynamical behavior [Eckmann et al., 1987; Webber Jr 



and Zbilut, 1994], e.g. diagonals (similar local time evolution of different 



parts of the trajectory) or horizontal and vertical black lines (state does not 
change for some time). 

Zbilut and Webber have developed the recurren c e qua ntification analysis 
(R QA) to quantify an RP [ ]Webber Jr and Zbiluj |1994 Zbilut and Webber 
Jr., 1992). They defined measures using the recurrence point density and 
diagonal structures in the recurrence plot, the recurrence rate RR (den- 
sity of recurrence points), the determinism DET (ratio of recurrence points 
forming diagonal structures to all recurrence points), the maximal length of 
diagonal structures L„iax (or their averaged length L), the Shannon entropy 
ENT of the distribution of the diagonal lengths and the trend TREND (pal- 
ing in the RP). The computation of these measures in shifted windows along 
the main diagonal of the RP enables one to find characteristic excursions of 
the trajectory in the phase space of the considered systems. 

Trulla et al. have applied these measures in order to find transitions in 
dynamical systems [ fTrulla et al. , 1996( ]. They have showed, that the RQA 
is able to find transitions between chaos and order (periodical states). But 
they could not find the chaos-chaos transitions. 



3 Laminarity and Trapping Time 

We have recently introduced two additional measures w hich are based on 
the vertical structures in the RP [ Marwan et al. , 2002p . We define these 



3 



measures analogous to the definition of DET and L, but we consider the 
distribution P{v) of the length of the vertical structures in the RP. 

First, the laminarity LAM 

LAM:=^^^^, (2) 

is the ratio of recurrence points forming vertical structures to all recurrence 
points and represents the probability of occurrence of laminar states in the 
system, but it does not describe the length of these laminar phases. It will 
decrease if the RP consists of more single recurrence points than vertical 
structures. 

Next, the trapping time TT 

rr^rr ._ T,v = 2 ^'-P(^) (o^ 
Ev=2 Pi'") 

is the averaged length of the vertical structures. The measure TT contains 
information about the amount and the length of the laminar phases. 

[Figure 1 about here.] 

The difference between these measures and the traditional RQA measures 



is their ability to find transitions between chaos and chaos [ Marwan et al 



2002). For example, such transitions can be found in the logistic map 
Xn+i = axn{l — Xn) with increasing control parameter a e [0,4] and .t„ g 
[0, 1] c TZ. For such trajectories x{a) which contain laminar states (e. g. a = 
3.678,3.791,3.927), LAM and TT show pronounced maxima (Fig. U). The 
application of these measures to heart rate variability data, has shown, 
that they are able to detect and quantify laminar phases before a life- 
threatening cardiac arrhythmia and, thus, to enable a prediction of such 



an event [ |Marwan et al.| , P002[ ]. These findings can be of importance for the 
therapy of malignant cardiac arrhythmias. 

In the next section we will apply this extended RQA to physiological data. 



4 Event Related Potentials 
4.1 The Oddball experiment 

As mentioned in the Introduction, the Oddball experiment studies brain 
potentials during a stimulus presentation. 

The measurement of the EEG was done with 31 electrodes/ channels (Tab. ||). 
The first 25 electrodes were localized as shown in Figure |2|; the others were 
reference electrodes. The sample interval for the measurements was 4 ms. 

[Figure 2 about here.] 
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[Table 1 about here.] 



Probands were seated in a dimly lit room in front of a monitor and were in- 
structed to count tones of high pitch. Each subject was tested in nine blocks. 
The blocks varied in the probability of occurrence of the higher tones from 
10 to 90 %. Each block contained at least 30 target tones. Response was 
given in a three alternative choice (using cursor keys of the keyboard). Dur- 
ing the test, the EEG was recorded. The stimuli were computer-generated 
beeps of 100 ms length. Tones were either high (1400 Hz) or low (1000 Hz). 
They were presented with an interstimulus interval of 1000 ms. 

After computing event-related voltage averages for the experimental ma- 
nipulations (10% up to 90% target probability) one can observe a P300 
ERF component whose amplitude is anti-correlated to the probability of 
the stimuli (surprise ERP, Fig. ^). 



[Figure 3 about here.] 



The P300 component of the ERP was the first potential discovered to vary 
in dependence on subject-internal factors like attention and expectation in- 
stead of physical characteristics [ Sutton et al. , 1965 1. The amplitude of the 



P300 component is highly sensitive to novelty of an event and its relevance. 
So this component is assumed to reflect the updating of the environmental 
model of the information processing system [context updating, Donchin , 



1981; Donchin and Coles, 198S]. 



4.2 Data analysis 

Our focus will be directed to the ERP data of two extreme event proba- 
bilities. Henceforth, the time (measured in ms) is denoted as t, the trial 
number as i and the electrode as e (the allocation of the electrode numbers 
with their notion, see Fig. 

The first set ERP90 contains 40 trials of ERP data for an event frequency of 
90 % and the second set ERPIO contains 31 trials for an event frequency of 
10 %. Both data sets can be rather well discriminated in the NlOO and P300 
components by the average over all trials (Fig. As expected, both com- 
ponents have increased for lower event probabilities (ERPIO). The maxima 
of the P300 are located around the central and central-parietal electrodes. 
However, the single trials do not obtain such a clear result. The P300 com- 
ponent is only well pronounced in 15 trials. When the single trials are 
observed, then extreme values can also occur in the ERP90 data and van- 
ish in the ERPIO data (Fig. ^). We applied also a statistical variance-based 
T-test to the single trial ERP data. However, this method could also not 
clearly distinguish the single trials. 



[Figure 4 about here.] 



The recurrence quantification analysis (RQA) is based on the structures 
obtained by recurrence plots (RPs). The RPs were firstly computed for the 
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means of ERP90 and ERPIO over all trials and then for the single trials. 
This was done with the embedding parameters m = 3, r = 3 and e = 10 % 
(fixed amount of nearest neighbours). The embedding parameters were es- 
timated by using the standard methods false nearest neighbours (dimen- 
sion) and mutual information (delay) [Kantz and Schreiber, 1997 1. Due 
to the NlOO and the P300 components in the data, the RPs show varying 
structures changing in time (Fig. |5|). Diagonal structures and clusters of 
black points occur. The nonstationarity of the data around the NlOO and 
P300 causes extended white bands along these times in the RPs. However, 
the clustered black points around 300 ms occur in almost all RPs of the 
ERPIO data set. 



[Figure 5 about here.] 



The RQA was computed from the RPs of ERP90 and ERPIO for the single 
trials, in sliding windows over the RPs (which have the dimension m = 3) 
with a length of 240 ms and with a shifting step of 8 ms. This window 
length corresponds with a data length of 60 values. 

The mean of all RQA variables of ERPIO reveal typical structures in the 
data (Fig. ^, right column). They indicate the transitions corresponding 
to the NlOO and P300 components around the central electrodes. The RQA 
variables for the ERP90 do not reveal these transitions (Fig. |^ left column). 
The onset of the increasing of the parameter is about 120 ms before the 
event. This is due to the windowed analysis of the RPs (240 ms windows). 
We have chosen the middle of the RP window for the time, what results in 
a 120 ms earlier onset of the RQA variables. 



[Figure 6 about here.] 



The four RQA variables are quite different, especially in their amplitude. 
For ERPIO, LAM and TT are the best pronounced parameter and have two 
distinct maxima at some electrodes; DET and L reveal these maxima at 
these electrodes too, but are lesser pronounced (Fig.|^. These maxima occur 
at the transition around 100 ms and 300 ms after the event and occur at the 
electrodes F3, F4, FZ, C3, FCZ, PZ, POZ and P03. Differences between the 
various transitions found by these measures also occur in time and brain 
locations (electrodes). But, the study was not detailed enough in order to 
give reliable results. 

The analysis of the single trials achieves similar results (Figs, [t] and ^ show 
the results for selected trials). The LAM clearly found the NlOO and P300 
components for ERPIO in 26 trials (of 31), but not for the ERP90 trials. The 
other measures have lesser maxima and, thus, are not suitable for such 
recognition. 

This result indicates that our introduced measures of complexity (espe- 
cially LAM) are able to recognize transitions in brain potentials, which are 
caused by e. g. stimulative events. These transitions can be found in the 
single trials, which is an improvement to the classical method of averaging 
all observations. 
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[Figure 7 about here.] 
[Figure 8 about here.] 

5 Summary 

We have apphed an extended recurrence quantification analysis (RQA) to 
physiological event related potential data (ERP). The classical RQA con- 
sists of measures which are mainly based on diagonal structures in the 
recurrence plots (RPs), e.g. the determinism (DET), which is the ratio of 
recurrence points located on connected diagonal structures in the RP, and 
the averaged diagonal line length (L). We have extended the RQA with 
two recently introduced measures, the laminarity (LAM) and the trapping 
time (TT). These measures are analogously defined as DET and L, but pro- 
vided by the vertical structures in recurrence plots. Whereas the classical 
RQA enables the identification of period-chaos transitions, the new mea- 
sures make the identification of chaos-chaos transitions and laminar states 
possible. 

The classical method to study ERP data is to average them over many trials. 
Our aim was to study the single trials in order to find transitions in the 
data. 

The application of the extended RQA to ERP data has discriminated the 
single trials with a distinct P300 component due to a high surprise moment 
(less frequent events) against such trials with a low surprise moment (high 
frequent events). Considering the raw ERPIO data, the P300 component 
can only be found in the half of all trials. Also a statistical variance test 
fails to distinguish cleary the trials. The LAM is the most pronounced pa- 
rameter in this analysis. It measures the ratio of recurrence points located 
on connected vertical structures in an RP. This structures correspond with 
laminarity within the underlying process. In the ERP data, the LAM re- 
veals transitions from less laminar states to higher laminar states after 
the occurrence of the event and a transition from higher laminar states 
to less laminar states after about 400 ms. These transitions occur around 
bounded brain areas (parietal to frontal along the central axis). The com- 
parable measures DET/ LAM and L/ TT are quite different in their am- 
plitude. There should also be differences in time and brain location of the 
found transitions. 

These results show that the measures based on vertical RP structures make 
the identification of transitions possible, which are not found by the clas- 
sical RQA measures. They indicate transitions in the brain processes into 
laminar states due to the surprising moment of observed events. 

A future work will be concerned with the development of a statistical eval- 
uation of these results. Furthermore, this investigation has to be extended 
to ERP data gained from other frequent events and a detailed study of the 
comparable measures DET/ LAM and L/ TT should give hints about the 
different transitions in the brain processes. 
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Figure 1: Laminarity (B) and trapping time (C) of time series gained from 
the logistic map for various control parameters (A). These measures reveal 
laminar and intermittent states. The vertical dotted lines show a choosing 
of points of band merging and laminar behaviour (a = 3.678, 3.727, 3.752, 
3.791, 3.877, 3.927). The length of the data were N = 1000 and the embed- 
ding parameters were = 1, t = 1 and £ = 0.1. 



12 




Figure 2: Localization of the electrodes on the head. 
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Figure 3: Mean event related potentials for event frequencies of 90 % (left, 
40 trials) and 10% (right, 31 trials). The NlOO and P300 components are 
well pronounced for the frequencies of 10 %. The lower plots show the ERP 
of selected electrodes. The reference of the electrode numbers is given in 
Table |. 
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Figure 4: Event related potentials for selected trials of the event frequen- 
cies of 90 % (left) and 10 % (right). Both, ERPIO and ERP90 of single trials 
can be strongly or weakly pronounced, respectively, which makes their dis- 
crimination difficult. The reference of the electrode numbers are given in 
Table |. 



15 



ERP90, trial 1 5 (CPZ) ERP1 0, trial 1 9 (CPZ) 




-200 200 400 600 800 -200 200 400 600 800 

Time [ms] Time [ms] 



Figure 5: ERP data for event frequencies of 90 % (upper left) and 10 % (up- 
per right), and their corresponding recurrence plots (lower plots). For the 
lower event frequency, more cluster of recurrence points occur at 100ms and 
300ms. 
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Figure 6: Averaged RQA measures for the ERP data of both event frequen- 
cies (averaged over all trials). Whereas the measures do not reveal any 
transitions in the ERP90 data, they clearly recognize the transitions for 
the ERPIO data. 



17 



ERP90, trial 15 (CPZ) 



ERP10, trial 19 (CPZ) 
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Figure 7: RQA measures for selected single trials and the central-parietal 
electrode (black). The trial-averaged RQA measures for the same electrode 
is shown in blue (the light blue band marks the 95 % significance interval). 
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Figure 8: RQA measures for the same 
electrodes. 
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Table 1: Notation of the electrodes and their numbering as it is used in the 
figures (electrodes 26-31 are reference electrodes). 



# 


Electrode 


# 


Electrode 


1 


F7 


14 


T8 


2 


FC5 


15 


P7 


3 


F3 


16 


PZ 


4 


FZ 


17 


P3 


5 


F4 


18 


CZ 


6 


FC6 


19 


P4 


7 


F8 


20 


P8 


8 


T7 


21 


oz 


9 


CPS 


22 


POZ 


10 


C3 


23 


P03 


11 


FCZ 


24 


CPZ 


12 


C4 


25 


P04 


13 


CP6 
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